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Details are presented of an efficient formalism for calculating transmission and reflection matrices 
from first principles in layered materials. Within the framework of spin density functional theory 
and using tight-binding muffin-tin orbitals, scattering matrices are determined by matching the 
wave-functions at the boundaries between leads which support well-defined scattering states and 
the scattering region. The calculation scales linearly with the number of principal layers N in the 
scattering region and as the cube of the number of atoms H in the lateral supercell. For metallic 
systems for which the required Brillouin zone sampling decreases as H increases, the final scaling 
goes as H 2 N. In practice, the efficient basis set allows scattering regions for which H 2 N ~ I0 6 to 
be handled. The method is illustrated for Co/Cu multilayers and single interfaces using large lateral 
supercells (up to 20 x 20) to model interface disorder. Because the scattering states are explicitly 
found, "channel decomposition" of the interface scattering for clean and disordered interfaces can 
be performed. 
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I. INTRODUCTION 

One of the most important driving forces in condensed 
matter physics in the last thirty years has been the con- 
trolled growth of layered structures so thin that inter- 
face effects dominate bulk properties and quantum size 
effects can be observed. In doped semiconductors, the 
large Fermi wavelength of mobile charge carriers made it 
possible to observe finite size effects for layer thicknesses 
on a micron scale. Much thinner layers must be used in 
order to make such observations in metals because Fermi 
wavelengths are typically of the order of an interatomic 
spacing. Nevertheless, following rapidly on the heels of 
a number of important discoveries in semiconductor het- 
erostructures, interface-dominated effects such as inter- 
face magnetic anisotropy, oscillatory exchange coupling 
and giant magnetoresistance (GMR) were found in arti- 
ficially layered transition metal materials. Reflecting the 
shorter Fermi wavelength, the characteristic length scale 
is of order of nanometers. 

Our main purpose in this paper is to give de- 
tails of a scheme we have developed which is 
suitable for studying mesoscopic transport in inho- 
mogeneous, mainly layered, transition metal mag- 
netic materials. In the context of a large num- 
ber of schemes designed to study transport either 
from first-principleo 1 ! 2 ! 3 ! 4 ! 5 ! 6 ! 7 ! 8 ! 9 ! 10 ! 11 ! 12 ! 13 ! 14 ! 15 ! 16 ! 17 . 18 or 
based upon electronic structures calculated from first- 
principlesJ^SSiSiiS^rfSiSi we will require our computa- 
tional scheme to be (i) physically transparent, (ii) first- 
principles, requiring no free parameters, (iii) capable of 
handling complex electronic structures characteristic of 
transition metal elements and (iv) very efficient in order 
to be able to handle lateral supercells to study layered 



systems with different lattice parameters and to model 
disorder very flexibly. A tight-binding (TB) muffin-tin- 
orbital (MTO) implementation of the Landauer-Buttiker 
formulation of transport theory within the local-spin- 
density approximation (LSDA) of density-functional- 
theory (DFT) will satisfy these requirements. 

Because wave transport through interfaces is naturally 
described in terms of transmission and reflection, the 
Landauer-Buttiker (LB) transmission matrix formulation 
of electron transport gained rapid acceptance as a pow- 
erful tool in the field of mesoscopic physics once the 
controversies surrounding the circumstances under which 
different expressions should be used had been resolved. 25 
The two-terminal conductance of a piece of material is 
measured by attaching leads on either side, passing a 
current through these leads and measuring the potential 
drop across the scattering region. In the LB formulation 
of transport theory, the conductance G is expressed in 
terms of a transmission matrix t = t(Ep) 

G = ^Tr{ttt} (1) 

where the element i M „ is the probability amplitude that a 
state \v) in the left-hand lead incident on the scattering 
region from the left (see Fig. ^) is scattered into a state 
\fi) in the right-hand lead. The trace simply sums over all 
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incident and transmitted "channels" v and ^ and 4- is the 
fundamental unit of conductance. In much current work 
on first-principles transport the conductance is calculated 
directly from Green's functions expressed in some conve- 
nient localized orbital representation, 27 Explicit calcula- 
tion of the scattering states is avoided by making use of 
the invariance properties of a trace. Because we want to 
make contact with a large body of theoretical literature 28 
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FIG. 1: Sketch of the configuration used in the Landauer- 
Biittiker transport formulation to calculate the two terminal 
conductance. A (shaded) scattering region (S) is sandwiched 
by left- (£.) and right-hand (TV) leads which have translational 
symmetry and are partitioned into principal layers perpendic- 
ular to the transport direction. The scattering region contains 
N principal layers but the structure and chemical composition 
are in principle arbitrary. 



on mesoscopic physics and address a wider range of prob- 
lems in the field of spin-dependent transport, we will cal- 
culate the microscopic transmission and reflection matri- 
ces t and r. By using a real energy, we will avoid the 
problems encountered in distinguishing propagating and 
evanescent states when a small but finite imaginary part 
of the energy is used. The Landauer-Buttiker formalism 
satisfies our first requirement of physical transparency. 

In developing a scheme for studying transport in tran- 
sition metal multilayers, a fundamental difference be- 
tween semiconductors and transition metals must be rec- 
ognized. Transition metal atoms have two types of elec- 
trons with different orbital character. The s electrons are 
spatially quite extended and, in solids, form broad bands 
with low effective masses; they conduct easily. The d 
electrons are much more localized in space, form narrow 
bands with large effective masses and are responsible for 
the magnetism of transition metal elements. The "mag- 
netic" electrons, however, being itinerant do contribute 
to electrical transport. The appropriate framework for 
describing metallic magnetism, even for the late 3d tran- 
sition metal elements, is band theory^* An extremely 
successful framework exists for treating itinerant electron 
systems from first-principles and this is the Local Den- 
sity Approximation (LDA) of Density Functional Theory 
(DFT). For band magnetism, the appropriate extension 
to spin-polarized systems, the local spin-density approx- 
imation (LSDA) satisfies our second requirement of re- 
quiring no free parameters^ 

Oscillatory exchange coupling in layered magnetic 
structures was discussed by Bruno in terms of generalized 
reflection and transmission matrices^* which were cal- 
culated by Stile a 32 i 33 i 34 for realistic electronic structures 
using a schem e 35 ' 3 ^ ' based on linearized augmented plane 
waves (LAPWs). At an interface between a non-magnetic 
and a magnetic metal, the different electronic struc- 
tures of the majority and minority spin electrons in the 



magnetic material give rise to strongly spin-dependent 
rencctioni^i Schep used transmission and reflection ma- 
trices calculated from first-principles with an embed- 
ding surface Green's function method— to calculate spin- 
dependent interface resistances for specular Co/Cu inter- 
faces embedded in diffusive bulk material^ The resulting 
good agreement with experiment indicated that interface 
disorder is less important than the spin-dependent reflec- 
tion and transmission from a perfect interface. Calcula- 
tions of domain wall resistances as a function of the do- 
main wall thickness illustrated the usefulness of calculat- 
ing the full scattering matrix*!* 3 ^ However, the LAPW 
basis set used by Stiles and Schep was computation- 
ally too expensive to allow repeated lateral supercells 
to be used to model interfaces between materials with 
very different, incommensurate lattice parameters or to 
model disorder. This is true of all plane-wave based ba- 
sis sets which typically require of order 100 plane waves 
per atom in order to describe transition metal atom elec- 
tronic structures reasonably well. 

Muffin-tin orbitals (MTO) form a flexible, minimal ba- 
sis set leading to highly efficient computational schemes 
for solving the Kohn-Sham equations of DFTi 40 i 41 i 42 i 43 
For the close packed structures adopted by the magnetic 
materials Fe, Co, Ni and their alloys, a basis set of 9 func- 
tions (s, p, and d orbitals) per atom in combination with 
the atomic sphere approximation (ASA) for the poten- 
tial leads to errors in describing the electronic structure 
which are comparable to the absolute errors incurred by 
using the local density approximation. This should be 
compared to typically 100 basis functions per atom re- 
quired by the more accurate LAPW method. MTOs thus 
satisfy our third and fourth requirements of being able 
to treat complex electronic structures efficiently. 

The tight-binding linearized muffin tin orbital (TB- 
LMTO) surface Green's function (SGF) method has been 
developed to study the electronic structure of interfaces 
and other layered systems. When combined with the 
coherent-potential approximation (CPA), it allows the 
electronic structure, charge and spin densities of lay- 
ered materials with substitutional disorder to be cal- 
culated self-consistently very efficiently^* In this paper 
we describe how we have combined a method for cal- 
culating transmission and reflection matrices based on 
wave- function matching (WFM), in a form given by 
Ando*^ for an empirical tight-binding Hamiltonian, with 
a first-principles TB-MTO basis. 42 Applications of the 
method to a number of problems of current interest 
in spin-transport have already been given in a num- 
ber of short publications: to the calculation of spin- 
dependent interface resistances where interface disorder 
was modelled by means of large lateral supercells^ to the 
first principles calculation of the so-called mixing con- 
ductance parameter entering theories of current-induced 
magnetization reversal*^ and Gilbert-damping enhance- 
ment via spin-pumping;— to a generalized scattering for- 
mulation of the suppression of Andreev scattering at a 
ferromagnetic/superconducting interface^ to the prob- 
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lem of how spin-dependent interface resistances influence 
spin injection from a metallic ferromagnet into a III-V 
semiconductor. 49 ' 50 These examples amply demonstrate 
that the fourth requirement is well satisfied. 

In Sec.m we give technical details of the formalism and 
illustrate it in Sec. IIIII where we calculate the transmis- 
sion matrices for clean and disordered Co/Cu interfaces, 
document a number of convergence and accuracy issues 
and give a detailed "channel-decomposition" analysis of 
the scattering in the presence of disorder. A comparison 
with other methods is made in Sec. IIVI 

II. THEORY 

Central to the wave-function matching method for cal- 
culating the transmission and reflection matrices is the 
equation of motion (EoM) for electrons with energy e, 
relating the vectors of coefficients Cj for layers I — 1, I, 
and I + 1: 

Hi,i-xCi-x + (H- e)i,iCi + n IJ+1 C I+1 = 0. (2) 

Here, Cj = Cn describes the wavefunction amplitude 
in terms of some localized orbital basis \i) of dimension 
M where i labels the atomic orbital and atom site. [For 
the muffin-tin orbitals to be outlined in Sec. Ill Al i will 
be a combined index Rim, where I and m are the az- 
imuthal and magnetic quantum numbers, respectively, of 
the MTO defined for an atomic-spheres-approximation 
(ASA) potential on the site R.] The EoM does not re- 
strict us to only considering nearest neighbour interac- 
tions since atoms can always be grouped into layers de- 
fined as to be so thick that the interactions between layers 
I and I ± 2 are negligible (see Fig. Such layers are 
called principal layers. Their thickness depends on the 
range of the interactions which in turn partly depends 
on the spatial extent of the orbital basis. It will be min- 
imized by using the highly localized tight-binding MTO 
representation. 

Consider the situation sketched in Fig. ^ where the 
scattering region S is contacted with left (£) and right 
(1Z) leads which have perfect lattice periodicity and sup- 
port well-defined scattering states. We assume that 
the ground state charge and spin densities and the cor- 
responding Kohn-Sham independent electron potential 
have already been calculated self-consistently. The cal- 
culation of the scattering matrix can now be split into 
two distinct parts. In the first stage, to be discussed in 
Sec. Ill Bl the eigenmodes of the leads (= Co for the /i- 
th mode), of which there are 2M, are calculated using an 
EoM appropriate to MTOs and making use of the lattice 
periodicity. By calculating their k vectors (which are in 
general complex) and velocities t>k, the eigenstates can be 
classified as being either left-going u /J (— ) or right-going 
u M (+). They form a basis in which to expand any left- 
and right-going waves and have the convenient property 
that their transformation under a lattice translation in 
the leads is easily calculated using Bloch's theorem (with 



k complex). We use the small Roman letters i,j to label 
the non-orthogonal basis and the small Greek letters /j, v 
to label the lead eigenmodes. 

In the second stage discussed in Sec. Ill CI a scatter- 
ing region S is introduced in the layers 1 < I < N which 
mixes left- and right-going lead eigenmodes. The scatter- 
ing region can be a single interface, a complex multilayer 
or a tunnelling junction, and the scattering can be in- 
troduced by disorder or simply by discontinuities in the 
electronic structure at interfaces. The v — ► fi element of 
the reflection matrix, r M „ , is defined in terms of the ratio 
of the amplitudes of left-going and right-going solutions 
in the left lead (in layer for example) projected onto 
the v th right-going and /I th left-going propagating states 
(k vector real) renormalized with the velocities so as to 
have unit flux. The scattering problem is solved by direct 
numerical inversion of a matrix with the leads included 
as a boundary condition so as to make finite the matrix 
which has to be inverted. 



A. Muffin Tin Orbitals and the KKR equation 

Muffin-tin orbitala 40 i 41 i 42 i 43 (MTO) are defined for po- 
tentials in which space is divided into non-overlapping 
atom-centred "muffin-tin" spheres inside which the po- 
tential is spherically symmetric and the remaining "in- 
terstitial" region where the potential is taken to be con- 
stant. The atomic spheres approximation (ASA) is ob- 
tained (i) by taking the kinetic energy in the intersti- 
tial region to be zero and (ii) by expanding the muffin- 
tin spheres so that they fill all space whereby the vol- 
ume of the interstitial region vanishes; for monoatomic 
solids such spheres are called atomic Wigner-Seitz (WS) 
spheres. Inside a WS (or MT) sphere at R, the solution of 
the radial Schrodinger equation regular at R, (j> R i(e,r R ) 
can be determined numerically for energy e and angular 
momentumum I resulting in the partial wave 

<t>Rlm{£,rn) = <pRL{s,r R ) = (j) R i{e, r R )Yi m (r R ) (3) 

where r R = r — R and r R = |r — R|. A continuous and 
differentiable orbital is constructed by attaching to the 
partial wave at the sphere boundary r R = s R a "tail" 
consisting of an appropriate linear combination of the 
solutions of the Laplace equation, 

J RL (r R )^ {r R /u) l [2{2l + l)]^Y L {r R ) (4) 

and 

K RL (r R ) = (r R M- l - l Y L (r R ), (5) 

which are respectively, regular at R and at infinity. u> is 
the average WS radius if the structure contains different 
atoms. In terms of the logarithmic derivative of <fii(e,r) 
at r = s 
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(4>'i(s, s) is the radial derivative), the radial solutions are 
matched if for r ^ s, 



l-Di 



Di-l 



21- 



)J?(r) 
(7) 



where we drop the explicit R-dependence when this does 
not give rise to ambiguity, or in terms of the potential 
function 



P°( £ ) = 2(2Z + 1)(^) 



Di{e) + l + l 
Di{e)-l 



and normalization Nf{e) = (f 



N?{e)4>i{e,r)=K*{r)-P?{s)J?{r) 



(8) 



(9) 



By subtracting from the partial wave, both inside and 
outside the MT sphere, the J rl {y r ) component which is 
irregular at infinity, a function is formed which is contin- 
uous, diffcrcntiable and regular in all space, an energy- 
dependent muffin tin orbital x°rl( £ ' v r) : 



X° RL (e,r R ) = N° Rl (e)cl> m (e,r) + P? 
= K° L (r R ) 



(10) 

tr Ss Sr 

(11) 



The tail K rl (yr) has the desirable property that closed 
forms exist for expanding it around a different site R' in 
terms of the regular solutions J R , L ,(r R >), 

K rl {vr) = -£ Jr'L'MSr, l , iRL (12) 



The expansion coefficients S R , L , RL form a so-called 
canonical structure constant matrix: they do not de- 
pend on the lattice constant, on the MT (or AS) po- 
tentials or on energy. Because of the augmentation with 
J rl (tr), the resulting MTO is no longer a solution of 
the Schrodinger equation (SE) inside its own sphere R. 
When, however, a solution of the SE is sought in the 
form of a linear combination of MTOs centred on differ- 
ent sites, 



*( £ ,r)=^ X O i ( £ ,r iJ )C° 



RL 



(13) 



R.L 



then the partial wave solution is recovered if the aug- 
menting term Jrl{yr) 011 s ^ e R ^ s cancelled by the tails 
of MTOs centred on all other sites R' ^ R, expanded 
about R. The condition for this to occur is the "tail- 
cancellation" condition: 



R'.L 



[P RL (e)6 R R,5 LL , - S RL>RIL ,] C% L , = 0. (14) 



All of the information about the structural geometry of 
the system under investigation is contained in the struc- 
ture constant matrix S RL R , L , while all of the information 
about the atomic species on site R needed to calculate 
the electronic structure (eigenvalues and eigenvectors) is 
contained in the potential functions Prl(£)- These are 
determined by solving the radial Schrodinger equation for 
the corresponding spherically symmetrical atomic sphere 
potential for energy e and angular momentum I. 

A disadvantage of these "conventional" MTOs is their 
infinite range. However, there is a remarkably simple gen- 
eralization of the MTOs which allows their range to be 
modified by introducing a set of "screening" constants 
am (not to be confused with the lead eigenmode in- 
dex) while the "tail-cancellation" condition remains es- 
sentially unchanged: 

E [PrlWrvSlv ~ Srl,wl>] C a R, L , = 0. (15) 

R!,V 



P a {e) is a diagonal matrix related to P°(e) by 

P°(e) 



P a (e) = P°(e) + P°(£)aP a (e) 



1 -aP°(s 



; (16) 



and 



S a = S° + S°aS a = S° (1 - aS°) \ (17) 



For any set of ccri , the energy-dependent MTOs with the 
normalization 



£~^)K^| 2 = i, 



R,L 



(18) 



form a complete set for the MT (AS) potential used 
in their construction. Here, P denotes an energy 
derivative and (|18f) follows from the relation N a (e) = 
[(u)/2)P a (e)] 1 ' 2 . Sets of parameters cxri have been found 
for which the "screened" structure constants S RL R , L , 
have very short range, decaying exponentially with the 
interatomic separation^ The set of parameters, (3 R i, 
which yields the shortest range MTOs is called the "tight- 
binding" (TB) representation^ For close-packed struc- 
tures, the range of S RL R , L , is in practice limited to first- 
and second-nearest neighbours. This TB set with a = (3 
is what we will use from now, unless stated otherwise, 
since it will allow us to define principal layers with a 
minimal thickness. 

For the determination of energy bands e(k), the tail- 
cancellation or KKR equations are inconvenient because 
the energy-dependence of the potential function makes 
it necessary to solve (|14|) or <|15[) by searching for the 
roots of a determinant, which is very time consuming. 
Much more efficient methods have been developed based 
on energy-independent MTOs. However, to study trans- 
port we only need to know P@(e) for a fixed energy, usu- 
ally the Fermi energy. We assume that the Kohn-Sham 
equations have already been solved self-consistently (us- 
ing for example a linearized method) so we have the po- 
tentials from which to calculate the potential functions. 
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Although (|15fl can be brought into Hamiltonian form by 
linearizing the energy dependent potential function (see 
Appendix , we will work directly with the more exact 
KKR equation. 

B. Eigenmodes of the leads 

We will assume that there exists two-dimensional 
translational symmetry in the plane perpendicular to 
the transport direction so that states can be charac- 
terized by a lateral wave vector k|| in the correspond- 
ing two-dimensional Brillouin zone. The screened KKR 
equation 4 -^ in the mixed representation of ky and real 
space layer index I (see Fig. is 

-Sj^iCi-i + (Pu(e) - S^) d - sf t » J+1 C/+i = 0, 

(19) 

where Cj = Cn = C IBlm is a (l max + l) 2 H = M di- 
mensional vector describing the amplitudes of the 7-th 
layer with H sites and (^ max + l) 2 orbitals per site. Pu 
and Si j are M x M matrices. Pjj is a diagonal matrix 
of potential functions characterizing the AS potentials of 
layer / and 

S Z= E ^(T)e ik "- T , (20) 
Te{Tj,j} 

where {Ti,j} denotes the set of vectors that connect one 
lattice site in the 7-th layer with all lattice sites in the 
J-th layer. 

By analogy with equation (|19f) is the equation of 
motion we will use to calculate the amplitudes of right- 
and left-going waves which determine the scattering ma- 
trix. We will solve it for a fixed value of e (usually ep), 
and some ky to find fe M (E, ky) the component of the Bloch 
vector in the transport direction. To keep the notation 
simple, explicit reference to the ky and s dependence 
will be omitted from now on. The formalism to be de- 
scribed in the following can be applied to any electronic 
structure code based on the KKR equation l|19|) , such as 
third-generation TB-LMTOi 51 i 52 i 53 

Let us first consider the Bloch states in the ideal 
lead. To obtain linearly independent solutions, we set 
C/ = A 7 Co, since in a periodic potential the wave func- 
tion should satisfy Bloch's theorem. The potential func- 
tion matrix is the same for all unit cells. The structure 
constant matrix depends only on the relative positions 
and, because that is how they are defined, there is only 
coupling between adjacent principal layers so the equa- 
tion of motion becomes 

/ So,l(P- S , ) -S^S lfi ^ C/ ) = A (g 7 \ 

1 < 21 ) 

The eigenvalue A M can be written in the form A M — 
exp(ik.^ ■ To) with To connecting equivalent sites in ad- 
jacent prinicpal layers. The wave vector k M can be de- 
composed into ky and a remainder which is in general 



not real, k M = (ky,k M - ky). Equation 1(21)1 has 2M 
eigenvalues and 2M eigenvectors, corresponding to M 
right-going and M left-going waves. By calculating the 
wavevectors and velocities (see Eq. (I37f) and Appendix lA"|) 
of the lead eigenmodes, the propagating and evanescent 
states can be identified and sorted into right-going or 
left-going modes. 

Letting ui (—),..., Um(—) denote the left-going solu- 
tions Co corresponding to eigenvalues Ai(— ),..., Am(— ) 
and ui(+), um{+) the right-going solutions corre- 
sponding to eigenvalues Ai (+),..., Am(+), the matrix 
Uifj,(±) is defined as 

U(±) = (ui(±)...u M (±)) (22) 

and the matrix A(±) as the diagonal matrix with ele- 
ments A /J (±). Following Ando, we next expand any left- 
or right-going wave, at I = for example, as 

C (±) = U(±)C(±). (23) 

Note that Co is a vector whose elements are labelled i 
while the elements of the vector C are labelled /i. 

F(±) = U(±)A ± U- 1 (±) (24) 

is the matrix of Bloch factors (including evanescent 
states) transformed onto the basis \i) and plays a cen- 
tral role in the following. Knowing it makes it possible 
to translate a state expressed in the basis |i) from layer 
J of the lead to layer / by 

C / (±)=F / - J (±)C J (±). (25) 

C. Scattering problem 

The scattering region <S, divided into N principal layers 
numbered 1 to N, is now inserted between the left and 
right leads. The resulting (scattering region + leads) 
problem is infinite dimensional in the real space MTO 
representation but by making use of their translational 
symmetry, the leads can be incorporated as boundary 
conditions and the scattering problem can be reduced to 
a finite problem whose dimension is determined by the 
size of the scattering region (number of sites x number 
of orbitals per site). 

We set about decoupling the scattering region from 
the leads, first on the left-hand side, then on the right. 
The amplitude in the 0-th layer is first separated into 
right- and left- going components Co = C (+) + C (— ). 
Because there is no additional scattering in the leads, 
the right- and left-going components can be translated 
to the left by one (principal layer) lattice spacing using 
the generalized Bloch factors (|25f) so the amplitude in 
layer —1 can be related to that in layer as 

C-i - F c 1 (+)C (+) + F c 1 (-)C (-) 
= [^ £ 1 (+)-F £ 1 (-)]Co(+) + ^ 1 (-)Co. (26) 
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allowing us to express C_i in terms of Co and Co(+) 
and so eliminate it from the equation of motion for the 
0-th layer 

— "So.-iC-i + (Po,o — So,o) Co — S'o.iCi = 0, (27) 

which becomes 

(-Po,o — So,o)Co — S'o^Ci 
= So.-x [Fc\+) - F c \-)] C (+). (28) 

Here C denotes the left lead and 5o,o = ^o.o + 
Sq,-\F2 {— )■ ~Sq,-!F^ {—) is the "embedding poten- 
tial" for the left lead and the net result is that the equa- 
tions of motion have been truncated at layer 0. 

I 



On the right-hand side of the scattering region, we are 
interested in the situation where only right-going waves 
can exist in the (N + l)-th layer, so 

C w+2 = Fn(+)C N+1 (+) (29) 

allowing us to eliminate Cn+2 from the EoM for Cn+i 

{Pn+i,n+i — Sn+i,n+i)Cn+i — Sn+i,nCn = 0, (30) 

where Sn+i,n+! = SV+i.tf+i + SW+i.JV+a-Fk (+) and 
— 6'jv+i ) 7v+2^7?,(+) is the embedding potential for the 
right lead. 



Making use of the lead boundary conditions, the tail cancellation condition for the scattering problem in real space 
is given by the set of inhomogeneous linear equations 



( {P- 5)o,o -Sb,i 
—Si,o (P — S)i,i —Si,2 

o -s 2A (p-s) 2 , 2 



\ 



(P-S) 



N,N 



/Co 
Ci 

Co 



Cat 

V Cjv+i / 



- (p-s) 



/Co 
Ci 

c 2 



Cjv 
V Cat+i 





— Sn.n+i 
Sn+i,n (P — S)n+i,n+i ) 

\ / So,-! [F c \+) - Fc\-)] C (+) \ 





(31) 



which can be solved in terms of g = ( P— S 



/C N 

Ci 

c 2 



V Cat+i / 



/So,-! [F^(+)-F^(-)] C (+)\ 



V 



/ 



This treatment is very similar to the widely used surface Green function method^ The boundary conditions in 1|31|) 
are explicitly defined by considering the Bloch wave coming from the left-hand side while for conventional retarded 
or advanced Green functions the boundary conditions are specified by an infinitesimal imaginary part of the energy 
parameter e. 

We are now in a position where we can relate the outgoing wave amplitude in the right electrode to the incoming 
wave in the left electrode through the Green function by 



C«+i(+) - Cjv+i - 9n+!,oSo,-! [F c \+) - F^(-)] C (+). 



(32) 



Using the transformation between the eigenstates and the localized basis functions [7j a (±), we obtain the transmission 
and reflection matrix elements^ 



^) 1/2 {U^(+) 9N+ !, So,-! [F c 1 (+)-F^(-)] U C (+)} 



(33) 
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^ = (^) V2 {^R <5o,o^o,-i [F c 1 (+] FZ\-)] 1) Uc(+)}^ , (34) 

where and v label Bloch states and , t>^ are the components of the corresponding group velocities in the transport 
direction. Similarly, an incident wave from the right side is transmitted or reflected as 

/ \ 1/2 

C = i U c \-)9o,n+iS n +i^ + 2 [FnH - F n (+)} Uni-)}^ (35) 



/ \ 1/2 



The group velocities in l|33|) - (|36[) are determined using the expression 

v,{±) = ^ [u], (±)5ft +1 u M (±)A M - h.c. 



(36) 



(37) 



which is derived in Appendix A. Here, I and I + 1 denote neighbouring principal layers in either left or right lead, 
d = To • n is the distance between equivalent monolayers in adjacent principal layers and n is a unit vector in the 
transport direction. 



D. Disorder 

Interfaces between materials with different lattice 
parameters^ and disordered interfaces'^ can be mod- 
elled very flexibly using lateral supercells. This approach 
allows us to study the effect of various types of disor- 
der on transport properties, ranging from homogeneous 
interdiffusion (alloying) to islands, steps etc. The su- 
percell description of disorder becomes formally exact in 
the limit of infinitely large supercells. In practice, sat- 
isfactory convergence is achieved for supercells of quite 
moderate size (see Sec. IIII C|) . 



Leads 

The factor limiting the "size" of scattering problem 
which can be handled in practice is the rank of the blocks 
of the block-tridiagonal scattering matrix in l|19l) . which 
is proportional to the number of atoms in the supercell. 
If performed straightforwardly in the manner outlined 
in Sec. Ill Bl the solution of the lead equation (|21|l in- 
volves solving a non-Hermitian eigenvalue problem whose 
rank is twice as large. Unless use is made of the greater 
translational symmetry present in the leads, this can be- 
come the limiting step in the whole calculation. Doing 
so makes it possible to reduce the dimension of the lead 
state calculation to a size determined by the dimension 
of a primitive unit cell which is usually negligible. 

We consider an Hi x iJ 2 lateral supercell defined by 
the real-space lattice vectors 

Ai = Hx&i and A 2 = H 2 a 2 (38) 



where ai and a 2 are the lattice vectors describing the 
in-plane periodicity of a primitive unit cell (Fig. The 
cells contained within the supercell are generated by the 
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FIG. 2: Illustration of lateral supercells and corresponding 
2D interface Brillouin zones. Top panel: lattice vectors for a 
primitive unit cell containing a single atom (lhs) and a 4 x 4 
supercell (rhs). Bottom panel: a single k-point in the BZ (rhs) 
corresponding to the 4x4 real-space supercell is equivalent to 
4x4 k-points in the BZ (lhs) corresponding to the real-space 
primitive unit cell. 
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set of translations 
T 



T e T = {T|| = ftiai + h 2 n 2 ; 

< hi < H u < h 2 < H 2 } (39) 



where T — 1, . . . , Hi x H 2 is a convenient cell index. In 
reciprocal space the supercell Brillouin zone is defined by 
the reduced vectors 



Bi = h x /H x and B 2 = b 2 /H 2 



(40) 



where bi and h 2 are the reciprocal lattice vectors cor- 
responding to the real space primitive unit cell. As a 
result the Brillouin zone (BZ) is folded down, as shown 
schematically in Fig. [5] (bottom rhs), and the single k| 
point (§ is used to label supercell quantities) in the su- 
percell BZ corresponds to the set of Hi x H 2 kn points 
in the original unfolded BZ 



{k,| = k| + hxBx + h 2 B 2 
0<hi<Hi,0<h 2 < H 2 } (41) 



with JC — 1, . . . , Hi x H 2 . Solutions associated with dif- 
ferent k||^- in the primitive unit cell representaton be- 
come different "bands" at the single k| in the supercell 
representaton. 

The indices T and K, provide a natural means of de- 
scribing the supercell-related matrices U s (±) and F s (±) 
and their inverses in terms of (Hi x H 2 ) 2 sub-blocks 
with dimensions defined by the primitive unit cell. Thus 
L/f-^i) is the block containing the amplitudes of the 
modes associated with k^ in the T-th real-space cell. 

Solving the single unit cell problem for the set of k||- 
points belonging to IK (lhs of Fig. |2J) and using the Bloch 
symmetry of the eigenmodes, we get trivially 



e ik n>c- T ii 



(42) 



where U (k|| K ) is the matrix \2'3 of modes for a primitive 
unit cell for k|| K and the ± qualifier has been dropped 
for simplicity. Defining the matrix of phase factors 



( 



\ 



Iff J 



(43) 



with H = Hi x H 2 , and its inverse Y = X 1 , we can 
straightforwardly determine 



c/ s (kf)l^ = c/- 1 (k ll jy Kr 



(44) 



and 



i^U( k f) = E X ^(k|| K )^T 2 ( 45 ) 



The procedure outlined above for determining the ma- 
trices describing the lead modes scales linearly with the 



size of the supercell i.e., as (Hi x H 2 ) rather than as 
(Hi x H 2 ) 3 which is the scaling typical for matrix opera- 
tions. Another advantage is that it enables us to analyse 
the scattering. By keeping track of the relation between 
supercell "bands" and equivalent eigenmodes at different 
k|| ^ (Fig. |2| we can straightforwardly obtain from i|33|) - 
(|36|l ^(ky^. , k|| K2 ) and other scattering coefficients. In 
other words the "interband" specular scattering in the 
supercell picture translates, in the presence of disorder 
in the scattering region, into the "diffuse" scattering be- 
tween the kn vectors belonging to the K set. 



III. CALCULATIONS 

Even though the theoretical scheme outlined above 
contains no adjustable parameters, its practical imple- 
mentation does involve numerous approximations, some 
physical, others numerical, which need to be evaluated. 
At present, any workable scheme must be based upon 
an independent particle approximation. The results of 
a transport calculation will be limited by the extent to 
which the single particle electronic structures used are 
consistent with the corresponding Fermi surfaces deter- 
mined experimentally using methods such as de Haas- van 
Alphen measurements or the occupied and unoccupied 
electronic states close to the Fermi energy determined 
by for example, photoclcctron spectroscopy. 

In this section we examine how various approximations 
affect our end results. We begin with the calculation 
of the scattering states in bulk Cu and bulk Co (fill Aj l. 
These are then used to study specular scattering from an 
ideal ordered Cu/Co(lll) interface (fill B|) after which we 
describe how we model disordered interfaces (|III Cll and 
how the results can be analysed. 



A. Leads 

For a crystalline conductor with Bloch translational 
symmetry, each state at the Fermi energy can move un- 
hindered through the solid so that the transmission ma- 
trix is diagonal with \t^ v \ = 5^ w . In this ballistic regime, 
Q reduces to 



G°(n) = -Y^K 



,(k|| 



(46) 



and calculation of the so-called Sharvin conductance be- 
comes a matter of counting the number of modes (chan- 
nels) propagating in the transport direction n, denoted 
in as N a (n). To solve (jJIJ in practice, the orbital 
angular momentum expansions in (|12|) and (|13|l . which 
are in-principle infinite, must be truncated by introduc- 
ing some cutoff in I, denoted Z max - Usually, a value of 
( mal — 2 or 3 is used, corresponding to spd- or spdf- 
bases. 
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FIG. 3: 

10 15 n _1 m 

ity spin) plotted as a function of the normalized area element 
used in the Brillouin zone summation, A 2 ]<l\\/Abz = 1/Q 2 . 
Q, the number of intervals along the reciprocal lattice vec- 
tor is indicated at the top of the figure. Squares represent 
the series (Q = 20, 40, 80, 160, 320) least-squares fitted by the 
dashed line; diamonds, the series (Q = 22, 44, 88, 176, 352) 
least-squares fitted by the dash-dotted line. The part of the 
curve for the Co minority spin case to the left of the vertical 
dotted line is shown on an expanded scale in the inset. An 
fee lattice constant of a = 3.6144 and spd basis were used 
together with von Barth-Hedin's exchange-correlation poten- 
tial. 



The k|| summation is carried out by sampling, on a reg- 
ular mesh, the 2D Brillouin zone (BZ) defined by the (lat- 
eral) translational periodicity perpendicular to n. The 
results of carrying out this BZ summation are shown in 
Fig.|3]where G a (h) is plotted as a function of A 2 1:\\/Abz, 
the normalized area element per ky -point for bulk fee Cu 
and for the majority and minority spins of bulk fee Co. 
When the 2D-BZ reciprocal lattice vectors are each di- 
vided into Q intervals, then A 2 k^/Asz = 1/Q 2 - It can 
be seen that the Sharvin conductance is converged to 
about 1% if 3600 = 60 x 60 points are used in the com- 
plete 2D-BZ and to about 0.2% for 102400 = 320 x 320 
sampling points. The worst case is for the minority spin 
of Co which has a complex multi-sheeted Fermi surface. 
To see if there are any simple underlying trends in the 
convergence, we repeatedly bisect the intervals used in 
the BZ summation starting with Q — 20 and Q — 22, 
shown in the figure as squares and diamonds, respectively 
and least-squares fitted with the dashed and dash-dotted 
lines. The convergence is fairly uniform but not very sys- 
tematic indicating that the summation is limited by fine 
structure in the integrand at the smallest length scale 
studied which can only be resolved by increasingly fine 
sampling. Thus there is nothing to be gained by develop- 
ing more sophisticated interpolation schemes and when 
we introduce disorder in Section IlII CI this will be even 
more so. However, in the following we will see that the 



level of convergence we can achieve with discrete sam- 
pling is quite adequate and not a limiting step in the 
whole procedure. 

The calculations shown in the figure were performed 
using an spd-basis, for an fee lattice constant a = 3.614 A 
corresponding to the experimental volume of bulk (fee) 
Cu and using the exchange-correlation potential calcu- 
lated and parameterized by von Barth and Hedin^ For 
convenience, and to avoid repetition, we will refer to 
this in the following as a "standard" configuration. The 
converged values are given (underlined) in Table [I] to- 
gether with values calculated using an fee lattice con- 
stant a — 3.549 A corresponding to the volume of bulk 
hep Com^ Because we shall be studying Cu/Co interfaces 
where the volume per atom is not known very precisely 
from experiment, we will want to estimate the variation 
that can be expected when different but equally reason- 
able lattice constants are used. The increase of 3.4% 
(from 0.558 to 0.577 x 10 15 Sl^m- 2 ) observed for Cu 
can be attributed to the increased areal density of Cu 
atoms, (3.614/3. 549) 2 corresponding to - 3.7%. The Ta- 
ble also contains the corresponding results obtained with 
an spdf-h&sis. To the numerical accuracy shown, there is 
no difference between the spd and spdf case for Cu. 

For Co majority spin states, there is a 4% decrease in 
the conductance on going from an spd to an spdf basis. 
For a lattice constant a = 3.614 A, the magnetization is 
1 .684 fiB I atom for an spd- and 1 .648 /i b /atom for an spdf 



basis corresponding, respectively, to n n 



5.342 and 



5.324 electrons in the majority spin bands. Since all five 
(nominal) majority-spin d bands are full there are 0.342 
and 0.324 electrons in the frcc-clcctron-like sp band. In 
a free electron picture the ratio of the projection of the 
spherical Fermi surfaces is (0.324/0.342) 2 / 3 = 0.96, thus 
explaining the observed numerical result. 

The Co majority-spin conductance scarcely changes 
with changing lattice constant, however. The origin of 
this behaviour lies in the volume dependence of the mag- 
netic moment. When the lattice constant is decreased, 
the d bands broaden and the magnetic moment decreases 
from 1.684 to 1.646 ^s/atom in the spd case with a corre- 
sponding decrease of the occupancy of the sp band from 
0.342 to 0.323 majority-spin electrons. The correspond- 
ing 4% decrease in conductance is almost perfectly com- 
pensated by the increased areal density of atoms so there 
is no net change. For the minority-spin conductance, the 
same factors play a role but now the d bands are only 
partly filled. This results in complex Fermi surfaces for 
which simple estimates cannot be made. In this case re- 
course must be made to full band structure calculations. 
We return to this in Sect. IIII Bl 

The calculations presented so far were carried out us- 
ing the exchange-correlation potential calculated and pa- 
rameterized by von Barth and HedinS This is only one 
of a number of potentials we could have used, none of 
which is clearly better than the others in describing the 
ground state properties of magnetic materials. To gauge 
the uncertainty arising from this arbitrary choice, a num- 
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G CT (111) 





o(A) 


basis 


n a 


present calc. 


Schep" 


Copper 


3.549 


spd 


5.5 


0.577(0.577,0.577) 


0.57 




3.549 


spdf 


5.5 


0.577(0.577) 






3.614 


spd 


5.5 


0.558(0.559) 


0.55 




3.614 


spdf 


•5.5 


0.558(0.558) 


0.55 


Cobalt 


3.549 


spd 


5.323 


0.469(0.459,0.467) 


0.45 


majority 


3.549 


svdf 


5.304 


0.449(0.440) 


0.43 




3.614 


spd 


5.342 


0.466(0.457) 


0.45 




3.614 


spdf 


5.324 


0.448(0.439) 




Cobalt 


3.549 


spd 


3.677 


1.082(1.081,1.082) 


1.10 


minority 


3.549 


spdf 


3.696 


1.120(1.125) 


1.13 




3.614 


spd 


3.658 


1.046(1.047) 


1.06 




3.614 


spdf 


3.676 


1.074(1.079) 





a Ref|33 

TABLE I: The Sharvin conductances per spin (in units of 
10 15 f2 _1 m~ 2 ) in the (111) direction for fee Cu and Co using 
the experimental volumes of Cu and Co. The underlined num- 
bers are the converged values discussed in relation to Fig. [3] 
Most of the results were obtained with von Barth-Hedin's 
exchange-correlation potential while the results in brackets 
are for Perdew-Zunger (PZ) and Vosko-Wilk-Nusair (VWN) 
parameterizations, respectively. Where a single number is 
given in brackets, it means that PZ and VWN potentials yield 
identical results to the accuracy given. The corresponding re- 
sults of Schep et al. are given in the last column. The number 
of electrons with spin a is given in the fourth column. 



ber of calculations were carried out using the potentials 
given by Perdew-Zunger— and Vosko-Wilk-Nusair— and 
the results are given in brackets in the Table. Using dif- 
ferent exchange-correlation potentials leads to variation 
in the conductances of the order of 1 or 2 percent. 

A different (but equivalent) approach was adopted by 
Schep et alm^ to the determination of the Sharvin con- 
ductances for the same systems using conventional first- 
principles LMTO-ASA bulk electronic band structures, 
i.e. using £i(k) rather than k^ie — £F,kii) as used here. 
He expressed the Sharvin conductance as a projection of 
the Fermi surface onto a plane perpendicular to the trans- 
port direction and calculated the areas using a suitably 
modified 3D-BZ integration scheme. His results are also 
given m Table [Q and are as consistent with our present 
values as can be expected when using two entirely differ- 
ent computer codes. 

In determining the conductance of the leads, the BZ 
summation does not present a problem. The uncertain- 
ties arising from small variations in the atomic volumes, 
from incompleteness of the basis and from the choice 
of LDA parameterization are of comparable size. The 
MTO-AS approximation can be systematically improved 
but only at substantial computational cost. Since there 
is currently no way to systematically improve upon the 
LDA we identify it and the lack of knowledge of the 
atomic structure as limiting factors in studying trans- 
port from first principles. Though the atomic structures 



could be determined theoretically by total energy mini- 
mization, the LDA again presents a barrier since it sys- 
tematically underestimates lattice constants of transition 
metals in particular of the 3d series. Gradient corrections 
sometimes yield improvements but unfortunately not sys- 
tematically so. We conclude that our knowledge of and 
ability to calculate from first principles Fermi surfaces 
for bulk magnetic materials such as Fe or Co does not at 
present justify using a more accurate but substantially 
more expensive computational scheme than the present 
one. 



B. Ordered Interfaces 

Cu and Co have slightly different atomic volumes. The 
equilibrium lattice constant of Cu is 3.614 A and of Co 
3.549 A, assuming an fee structure. Even in the absence 
of interface disorder, the lattice spacing will not be ho- 
mogeneous and will depend on the lattice constant of the 
substrate on which the sample was grown, on the global 
and local concentrations of Cu and Co, and on other de- 
tails of how the structure was prepared. In principle we 
could calculate all of this by energy minimization. How- 
ever, we judge that the additional effort needed is not 
justified by current experimental knowledge. Instead, we 
content ourselves with estimating the uncertainty which 
results from plausible variations in the (interface) struc- 
ture by considering two limiting cases and one interme- 
diate case. In each case an fee structure is assumed, with 
lattice constants corresponding to (i) the atomic volume 
of Cu, (ii) the atomic volume of Co, (iii) an intermediate 
case with arithmetic mean of Cu and Co atomic volumes. 



a(A) 


3.549 


3.582 


3.614 


Basis 


spdf 


spd 


spd 


spd 


mcu(bulk) 


0.000 


0.000 


0.000 


0.000 


mcu(int-4) 


0.001(1) 


0.001 


0.001 


0.001 


racu(int-3) 


-0.001(0) 


0.000 


0.000 


0.000 


mcu(int-2) 


-0.005(5) - 


-0.005(4, 5)- 


-0.005(4) 


-0.005 


mcu(int-l) 


0.002(4) 


0.004(6,4) 


0.003(4) 


0.001(2) 



mco(int+l) 
mco (int+2) 
mco (int+3) 
mco (int+4) 
mco(bulk) 



1.526(490) 1 
1.621(597) 1 
1.602(576) 1 
1.610(587) 1 
1.609(590) 1 



578(45, 73)1. 
656(35, 53)1. 
645(21,41)1. 
649(27, 45)1. 
646(22,42)1. 



605(573) 

673(53) 

662(39) 

665(45) 

667(45) 



1.636(01) 
1.690(70) 
1.680(59) 
1.683(62) 
1.684(62) 



G maj (lll) 
G min (lll) 



0.409(399) 0. 
0.378(379) 0. 



431(21,29)0. 
378(80, 79)0. 



433(22) 
371(73) 



0.434(24) 
0.364(67) 



TABLE II: Variation of the layer-resolved magnetic moments 

(in Bohr magnetons) for Cu/Co(lll) interfaces with basis 

set and lattice constant. These results were obtained with 

von Barth-Hedin's exchange-correlation potential while the 

results in brackets, where given, are for Perdew-Zunger and 

Vosko-Wilk-Nusair parameterizations, respectively. In the 

last two rows, the interface conductances are given in units of 
lis 
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Our starting point is a self-consistent TB-LMTO SGF 
calculation^ for the interface embedded between semi- 
infinite Cu and Co leads whose potentials and spin- 
densities were determined self-consistently in separate 
"bulk" calculations. The charge and spin-densities are 
allowed to vary in nc u layers of Cu and nco layers of Co 
bounding the interface. The results of these calculations 
for Cu/Co(lll) interfaces and the three different lattice 
constants detailed above are given in Table ITT1 for nrj u =4, 
7ic =4. In the Cu layers, only tiny moments are induced. 
Only four layers away from the interface on the Co side, 
the magnetic moments are seen to be very close to the 
bulk values. At the interface, where the d-bandwidth is 
reduced as a result of the lower coordination number, the 
moments are suppressed rather than enhanced. This oc- 
curs because the majority-spin d bands are full and their 
number cannot increase. The width of the free-electron 
like sp band is less sensitive to the change in coordina- 
tion and its exchange splitting also changes less. As a 
result, there is little change in the sp moment. When the 
c?-band width is reduced, there is conversion of minority- 
and majority-spin sp electrons, without loss of the sp mo- 
ment, to the minority-spin d band with loss of d moment. 
This picture is supported by the full calculations. 

Earlier we saw that an ~ 2% change in lattice con- 
stant changed the bulk magnetic moment of fee Co by 
2.3%. The effect of changing the basis, from spd to spdf, 
was similar. From Table ^ the interface moments are 
seen to behave in a comparable fashion. The magnetic 
moment of the interface Co atoms decreases by 3.7%, 
from 1.636 yus/atom for a — 3.614 A to 1.578 (j.b /atom 
for a = 3.549 A for an spd basis and decreases from 
1.578 /is/atom to 1.526 /is/atom for an spdf basis for 
a = 3.549 A, a change of 3.4%. Thus the sp to d m i n 
conversion is enhanced at the interface by the reduced 
(i-bandwidth. 

Once the interface potential has been obtained, the 
transmission matrix can be calculated and the BZ sum- 
mation carried out. The convergence of this summation, 
shown in Fig. 0] for a lattice constant of a = 3.614 A and 
an spd basis closely parallels that seen in Fig.[3]and there- 
fore the k-summation does not represent a limitation in 
practice. Converged conductances 



At,!/,k|i 



h 



KuW (47) 



are given in the last two rows of Table [n] Though we 
will not concern ourselves in this publication with the 
application of the formalism we have been developing to 
a detailed interpretation of experimental observations, it 
should be noted that even a modest spin-dependence of 
"bare" interface conductances (~ 20%) can lead to spin- 
dependent interface resistances differing by a factor of 
^3 — 5 once account is taken of the finiteness of the 
conductance of the perfect leads using a formula derived 
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FIG. 4: Interface conductance G CT (111) (in units of 
10 15 _1 m" 2 ) for an fee Cu/Co(lll) interface for major- 
ity and minority spins plotted as a function of the normal- 
ized area element used in the Brillouin zone summation, 
A 2 k||/yl_B2 = 1/Q 2 . Q, the number of intervals along the 
reciprocal lattice vector is indicated at the top of the fig- 
ure. Squares represent the series (Q = 20, 40, 80, 160, 320) 
least-squares fitted by the dashed line; diamonds, the series 
(Q = 22, 44, 88, 176, 352) least-squares fitted by the dash- 
dotted line. The part of the curve for the Co minority spin 
case to the left of the vertical dotted line is shown on an ex- 
panded scale in the inset. A "standard" configuration was 
usedi£2 



by Schep et ali 
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(48) 



where and Ng are the Sharvin conductances (see Eq. 
(|46|l ) of the materials A and B forming the interface, in 
units of e 2 /h. 

The majority-spin case can be readily understood in 
terms of the geometry of the Fermi surfaces of Cu and 
Co so we begin by discussing this simple case before ex- 
amining the more complex minority-spin channel. 



Clean Cu/Cu (111) Interface: Majority Spins 

In the absence of disorder, crystal momentum parallel 
to the interface is conserved. If, for a given value of kii , 
there is a propagating state in Cu incident on the inter- 
face but none in Co, then an electron in such a state is 
completely reflected at the interface. Conversely, kii 's for 
which there is a propagating state in Co but none in Cu 
also cannot contribute to the conductance. To determine 
the existence of such states, it is sufficient to inspect pro- 
jections of the Fermi surfaces of fee Cu and majority-spin 
Co onto a plane perpendicular to the transport direction 
n, shown in Fig.[5]for n — (111). The first feature to note 
in the figure (left-hand and middle panels) is that per ky 
there is only a single channel with positive group velocity 
so that the transmission matrix in (|47|) is a complex num- 
ber whose modulus squared is a transmission probability 
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FIG. 5: Top row, left-hand panel: Fermi surface (FS) of Cu; 
middle panel: majority-spin FS of Co; right-hand panel: Cu 
FS viewed along the (111) direction with a projection of the 
bulk fee Brillouin zone (BZ) onto a plane perpendicular to 
this direction and of the two dimensional BZ. Bottom row, 
left-hand and middle panels: projections onto a plane per- 
pendicular to the (111) direction of the Cu and majority-spin 
Co Fermi surfaces; right-hand panel: transmission probabil- 
ity for majority-spin states as a function of transverse crystal 
momentum, T(kii) for an fee Cu/Co(lll) interface. A "stan- 
dard" configuration was used^ 



with values between and 1. It is plotted in the right- 
hand panel and can be interpreted simply. Regions which 
are depicted blue correspond to k|| 's for which there are 
propagating states in Cu but none in Co. These states 
have transmission probability and are totally reflected. 
For values of k|| for which there are propagating states 
in both Cu and Co, the transmission probability is very 
close to one, depicted red. These states are essentially 
free electron-like states which have the same symmetry 
in both materials and see the interface effectively as a 
very low potential step. Close to the centre of the figure 
there is an annular region where there are propagating 
states in Co but none in Cu so they do not contribute to 
the conductance. Performing the sum in l|47|l . we arrive 
at an interface conductance of 0.434 x 10 15 Q _1 m~ 2 to 
be compared to the Sharvin conductances given in Ta- 
ble Q] for Cu and Co; for a = 3.614 A and an spd basis 
these are, respectively, 0.558 and 0.466 in the same units. 
The interface conductance of 0.434 is seen to be essen- 
tially the Sharvin conductance of the majority states of 
Co reduced because the states closest to the A-axis (cor- 
responding to the symmetry axis of the figures, the TL 
line in reciprocal space) do not contribute. The explana- 
tion of the 5% decrease found on going from an spd to 
an spdf basis, (0.431 to 0.409), parallels that given for 
the corresponding change in the Sharvin conductance of 
bulk Co (0.469 to 0.449 in Tabled). 



Clean Cu/Cu (111) Interface: Minority Spins 

The minority-spin case is considerably more complex 
because the Co minority-spin d bands are only partly 
filled, resulting in multiple sheets of Fermi surface. These 
sheets arc shown in Fig. Altogether with their projections 
onto a plane perpendicular to the (111) transport direc- 
tion. Compared to Fig. [5] one difference we immediately 
notice is that even single Fermi surface (FS) sheets are 
not single valued: for a given k|| there can be more than 
one mode with positive group velocity. The areas de- 
picted green in the projections of the FS sheets from the 
fourth and fifth bands are examples where this occurs. 

An electron incident on the interface from the Cu side, 
with transverse crystal momentum ky , is transmitted 
into a linear combination of all propagating states with 
the same k|| in Co: the transmission matrix ^^(ky) is 
in general not square but rectangular. The transmission 
probabilities T p;y (k||) are shown in the bottom row of 
Fig. Because there is only a single incident state for all 
k|| , the maximum transmission probability is one. Com- 
parison of the total minority-spin transmission probabil- 
ity Tcft(k||) (Fig. bottom right-hand panel) with the 
corresponding majority-spin quantity (right-hand panel 
of Fig.JSJ strikingly illustrates the spin-dependence of the 
interface scattering, much more so than the integrated 
quantities might have led us to expect; the interface con- 
ductances, 0.364 and 0.434 x lO^fTW 2 from TableHU 
differ by only ~ 20%. 

Three factors contribute to the large k|| -dependence 
of the transmission probability: first and foremost, the 
complexity of the Fermi surface of both materials but 
especially of the minority spin of Co; secondly and inex- 
tricably linked with the first because of the relationship 
hvis. = Vfce(k), the mismatch of the Fermi velocities of 
the states on either side of the interface. Thirdly, the or- 
bital character of the states /x and v which varies strongly 
over the Fermi surface and gives rise to large matrix ele- 
ment effects. 

The great complexity of transition metal Fermi sur- 
faces, clear from the figure and well-documented in stan- 
dard textbooks, is not amenable to simple analytical 
treatment and has more often than not been neglected 
in theoretical transport studies. Nevertheless, as illus- 
trated particularly well by the ballistic limitfi^ spin- 
dependent band structure effects have been shown to 
lead to magnetoresistance ratios comparable to what are 
observed experimentally in the current-perpendicular-to- 
plane (CPP) measuring configuration and cannot be sim- 
ply ignored in any quantitative discussion. Most at- 
tempts to take into account contributions of the d states 
to electronic transport do so by mapping the five d bands 
onto a single tight-binding or free-electron band with a 
large effective mass. 

Fermi surface topology alone cannot explain all aspects 
of the tranmission coefficients seen in Fig. For exam- 
ple, there are values of k|| , such as that labelled Y in the 
figure, for which propagating solutions exist on both sides 
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FIG. 6: Top row, lefthand panel: Fermi surface (FS) of fee Cu; middle panels: third, fourth and fifth FS sheets of minority-spin 
fee Co; righthand panel: projection of the bulk fee Brillouin zone (BZ) onto a plane perpendicular to the (111) direction and 
of the two dimensional BZ. Middle row: corresponding projections of individual FS sheets and (rhs) of Co total. The number 
of propagating states with positive velocity is colour-coded following the colour bar on the right. Bottom row: probability 
T M „(k||) for a minority-spin state on the single FS sheet of Cu (y = 1) to be transmitted through a Cu|Co(lll) interface into 
FS sheet \l of fee Co as a function of the transverse crystal momentum ky . The point Y is such that there are only propagating 
states in Cu and in the fourth FS sheet of Co. For the point Y' slightly further away from A and indicated by a small open 
square there is, in addition, a propagating state in the third FS sheet of Co. Results are for a "standard" configuration^- 



of the interface yet the transmission probability is zero. 
This can be understood as follows. At k|| = Y, the prop- 
agating states in Cu have {s,p y ,p z ,d yz ,d 3z 2_ r 2,d x 2_ y 2} 
character (assuming the choice of in-plane axes as illus- 
trated in the top righthand panel of Fig. |5J) and are even 
with respect to reflection in the plane formed by the 
y-axis and the transport direction perpendicular to the 
(111) plane which we choose to be the z-axis. For this k|| 
the only propagating state in Co is in the fourth band. It 
has {p x , d xy , d xz } character which is odd with respect to 
reflection in the yz plane. Consequently, the correspond- 
ing hopping matrix elements in the Hamiltonian (and in 
the Green's function) vanish and the transmission is zero. 

Along the k y axis the symmetry of the states in Cu 
and those in the fourth band of Co remain the same and 
the transmission is seen to vanish for all values of k y . 
However, at points further away from A, we encounter 
states in the third band of Co which have even character 
whose matrix elements do not vanish by symmetry and 
we see substantial transmission probabilities. Similarly, 
for points closer to A, there are states in the fifth band 
of Co with even character whose matrix elements also 



do not vanish and again the transmission probability is 
substantial. Because it is obtained by superposition of 
transmission probabilities from Cu into the third, fourth 
and fifth sheets of the Co FS, the end result, though it 
may appear very complicated, can be straightforwardly 
analysed in this manner k-point by k-point. 



Though the underlying lattice symmetry is only three- 
fold, the Fermi surface projections shown in Fig. have 
six-fold rotational symmetry about the line A because 
the bulk fee structure has inversion symmetry (and time- 
reversal symmetry). The interface breaks the inversion 
symmetry so T^fa^) has only threefold rotation symme- 
try for the individual FS sheets. However, in-plane in- 
version symmetry is recovered for the total transmission 
probability Tcti{— ki|) = 7£-R.(k||) which has full sixfold 
symmetry. This follows from the time-reversal symmetry 
and is proven in Appendix 151 
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FIG. 7: Interface conductance (in units of 10 15 Q 1 m 2 ) for 
a disordered Cu/Co (111) interface modelled as 2ML of 50%- 
50% alloy in a yH x \fH lateral supercell as a function of 
\fS. The results are given for different randomly generated 
configurations of disorder (15 for minority spin, 5 for majority 
spin). Results are for a "standard" configuration!— 

C. Interface Disorder 

Instructive though the study of perfect interfaces may 
be in gaining an understanding of the role electronic 
structure mismatch may play in determining giant mag- 
netoresistive effects, all measurements are made on de- 
vices which contain disorder, mostly in the diffusive 
regime. Because there is little information available from 
experiment about the nature of this disorder, it is very 
important to be able to model it in a flexible manner, 
introducing a minimum of free parameters. To model in- 
terfaces between materials with different lattice constants 
and disorder, we use the lateral supercells described in 
section Hi Dl Since this approach is formally only valid if 
sufficiently large supercells are used, we begin by study- 
ing how the interface conductance depends on the lateral 
supercell size. 

To perform fully self-consistent calculations for a num- 
ber of large lateral supercells and for different configura- 
tions of disorder would be prohibitively expensive. Fortu- 
nately, the coherent potential approximation (CPA) is a 
very efficient way of calculating charge and spin densities 
for a substitutional disordered A X B\- X alloy with an ex- 
pense comparable to that required for an ordered system 
with a minimal unit celliSi The output from such a calcu- 
lation are atomic sphere potentials for the two sites, va 
and vb- The layer CPA approximation generalizes this 
to allow the concentration to vary from one layer to the 
next 4^ 

Once va and vb have been calculated for some con- 
centration x, an H = H\ x iJ 2 lateral supercell is con- 
structed in which the potentials are distributed at ran- 
dom, maintaining the concentration for which they were 
self-consistently calculated. The conductances calculated 
for 4 < VH < 20 are shown in Fig. □ for a Cu/Cu(lll) 



interface in which the Cu and the Co layers forming the 
interface are totally mixed to give two layers of 50%-50% 
interface alloy. For each value of H, the results for a 
number of different randomly generated disorder config- 
urations are shown (20 for minority, 5 for minority spin). 
The sample to sample variation is largest for the minority 
spin case, ranging from ±5% for a modest 4x4 unit cell 
and decreasing to less than ±1% for a 20 x 20 unit cell. 
For \[H ~ 10, the spread in minority spin conductances 
is ~ 5% which is comparable to the typical uncertainty 
we associated with the LDA error, the uncertainty in lat- 
tice constants or the error incurred by using the ASA. 

Comparing now the conductances without and with 
disorder, we see that interface disorder has virtually 
no effect on the majority spin channel (0.434 versus 
0.432 x 10 15 f2 _1 m~ 2 ) which is a consequence of the 
great similarity of the Cu and Co majority spin potentials 
and electronic structures. However, in the minority-spin 
channel the effect (0.364 versus 0.31 x 10 15 Q -1 m -2 ) is 
much larger. As noted in the context of (|48|l . a relatively 
small change in the interface transmission can lead to 
a large change in the interface resistance when account 
is taken of the finite conductance of the leads. We will 
return to the consequences for the spin-dependent inter- 
face resistance after completing the study of the interface 
transmission on which it is based. 

When disorder is modelled in lateral supercells, the 
transmission probabilities can be classified as being spec- 
ular or diffuse depending upon whether or not transverse 
momentum is conserved^Si In the presence of interface 
disorder, the conductance per unit area can be expressed 
as 

G = G 8 + Gd 

2 2 

= 7rE T ^( k ih k n) + ir E T /-( k u>4) ( 49 ) 

k ll kj^kj, 

where k|| and k'| belong to the two dimensional Bril- 
louin zone for (1 x 1) translational periodicity and 
^(k^kji) = t (UI/ (k||,kj|)tt I/ (k||,kj|). The transmission 
matrix elements between two Bloch states with the same 
k|| are defined to be specular, those between scattering 
states with different k|| as being diffuse. In the absence of 
interface disorder, there is by definition only a specular 
component. 

Dependence of interface conductance on alloy concentration 

The results in Fig. were obtained for a structural 
model of the Co/Cu(lll) interface consisting of two 
monolayers (2ML) of 50%-50% alloy that was derived 
from X-rayS NMR&Gi and magnetic EXAFSS studies. 
Though the most plausible model there is at present, it 
contains large uncertainties. This makes it important to 
explore the consequences of varying the parameters defin- 
ing the model. To do so, we calculate the conductance 
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FIG. 8: Illustration of 3 different models of interface dis- 
order considered. Top (1ML): disorder is modelled using 
one monolayer (ML) of [Cui-^Coa;] alloy between Cu and 
Co leads, denoted as CufCui-xCo^Co. Middle (2ML): dis- 
order modelled in two MLs as Cu[Cui_ x Co x |Cu2,Coi_a;]Co. 
Bottom (4ML): starting from the 2 ML disorder case, 
1/3 of the concentration x of impurity atoms is trans- 
ferred to the next layer resulting in disorder in four MLs: 
CufCux.x Cosl |Cuj_ 2x Co2x |Cu2 £ Co 1 _ 2sl | Cu i!i Cox_ x ] Co. 
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FIG. 9: Interface conductance of a disordered Cu/Co (111) 
interface with disorder modelled in a 20 x 20 lateral supercell 
as a function of the alloy concentration x. Results are shown 
for the three different models described in Fig.|8]with disorder 
in 1, 2 or 4 MLs. Only a single disorder configuration was 
used and the size of the symbols corresponds to the spread 
in values found for this supercell size in Fig. |7| For 1ML, 
the total conductance is resolved into specular and diffuse 
components (see the legend). Results are for a "standard" 
configuration^ 



using 20x20 lateral supercells as a function of alloy con- 
centration for models in which the disorder is confined 
to one, two or four monolayers. The three models are 
defined in Fig. [SJ From the results shown in Fig. |5J it 
can be seen that the interface transmission for majority- 
spin electrons depends only very weakly on alloy concen- 
tration and its spatial distribution. The results for the 
1ML, 2ML and 4ML models cannot be distinguished on 
the scale of the figure. When the conductance is decom- 
posed using Ij49(l . the diffuse component is found to be 
very small. Therefore, only the results for the minority- 
spin case need be examined in any detail. 

We start by varying the alloy concentration over the 
full concentration range (0-100%) for a disordered mono- 
layer. The magnitude of resulting variation in transmis- 
sion is limited (~ 7%) but exceeds the statistical spread 
between different configurations of disorder, which ac- 
cording to Fig. [7| is less than ±1%. On adding Co to a 
layer of Cu, the transmission decreases, reaches a min- 
imum for 10% Co, then increases monotonically up to 
70-90% region where the transmission is higher than for 
a clean interface. 66 100% Co represents a clean interface 
again, so this limit must yield the same transmission as 
0% Co. 

The variation can be examined in terms of the specular 
and diffuse components defined in 149(1 . From Fig. EI it 
can be seen that, for the minority spin channel, the dif- 
fuse scattering by Co impurity atoms in Cu is stronger 
than that by Cu impurity atoms in Co. However, the 
specular scattering is also more strongly reduced by Cu 
in Co than by Co in Cu. The two effects largely cancel 
resulting in the observed undulatory total transmission 
as a function of the alloy concentration seen in the fig- 



ure. The diffuse scattering has a maximum close to a 
50%-50% alloy concentration where its contribution to 
the conductance is almost twice as large as from the 
specular scattering. While the conductance as such is 
scarcely effected, the strong diffuse scattering will play 
an important role in destroying the phase coherence of 
the electrons. Ultimately, this will be the physics under- 
lying the so-called "two current series resistor" (2CSR) 
modelM^ 

If the disorder extends over more than a monolayer, 
then modelling the interface as several layers of homo- 
geneous alloy is not obviously realistic. Instead, one 
might expect the layers closest to the interface to be 
most strongly mixed, the amount of mixing decreasing 
with the separation from the interface. A simple way to 
model this is to take two interface layers, one Cu and 
one Co, and to mix them in varying degrees. Denoting 
this Cu|Co interface as CulCui-xCo^lCuzCoi-JCo we 
consider < x < 0.5 i.e., the Cu concentration decreases 
monotonically from left to right. The calculated interface 
transmission is seen to essentially interpolate linearly the 
results obtained previously for the clean (x — 0) and dis- 
ordered (x = 0.5) cases. 

A slightly more elaborate model can be constructed 
from the 2ML model by distributing the x impurity 
atoms so that 2x/3 are in the interface layer while x/3 
are to be found further from the original interface, in the 
following layer. This results in the concentration profile 
Cu[Cui_*Co.e|Cui 2* C02* |Cu2x Co, 2« |CuaCoi_»lCo. 

L 3 3 1 i_ -3" "3" 1 ~3~ L ~—' 3 3 1 

x = corresponds to a completely ordered interface 
while the maximum value x can have so that the 
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concentration decreases from left to right monotonically 
is 75%. This relatively small redistribution of intermixed 
atoms is seen to reduce the transmission by 15% for 
x = 0.5. Even for very good metals, relatively opaque 
interfaces can result when the electronic structure on 
either side have different characters. In such situations, 
disorder can influence the transmission strongly even 
reducing the interface resistance very substantially^ A 
detailed analysis of the different contributions to the 
interface scattering in the 2ML and 4ML cases will be 
given in a separate publication. 



D. Analysis of Interface Disorder Scattering 

The scattering induced by two layers of 50%-50% alloy 
is illustrated in Fig. 1101 and Fig. ^] for the majority and 
minority spins, respectively, of a Cu/Co(lll) interface. 
Calculations were performed for the single kjj point, T, 
and a 20 x 20 lateral supercell equivalent to using a 1 x 1 
interface cell and k-space sampling with 20 x 20 points 
in the corresponding BZ. Disorder averaging was carried 
out using 5 (for majority spin) or 20 (for minority spin) 
disorder configurations generated randomly. 

Figs. HOf a') and llOf b) show the majority-spin Fermi 
surface projections of fee Cu and Co, respectively, ob- 
tained from "unfolding" the supercell calculation. The 
coarse 20 x 20 grid is seen to yield a good representa- 
tion of the detailed Fermi surface projections shown in 
Fig. T(k|| , kf, ) is shown in Fig. EJc) for k|| = Y 




FIG. 10: Fermi surface projections of majority-spin fee Cu 
(a) and Co (b) derived from a single k-point using a 20 x 20 
lateral supercell. The dark red point in the Cu Fermi surface 
projection corresponds to the point Y in the top righthand 
panel of Fig.HJ T(Y, kji) is shown in (c), and in (d) magnified 
by a factor 500 where the ballistic component T(Y, kn = Y) 
is indicated by a white point because its value goes off the 
scale. Results are for a "standard" configuration— and aver- 
aged over 5 different configurations of disorder. 



on the k y axis in Fig. Specular scattering dominates 
with T(k|| = Y, kf, = Y) = 0.93. The diffuse scatter- 
ing is so weak that nothing can be seen on a scale of T 
from to 1. To render it visible, a magnification by a 
factor 500 is needed, Fig. llOf dh The total diffuse scat- 
tering, T d (Y) = T(k„ = Y,k\ / 7) = 0.04 can 

be seen from the figure to be made up of contributions 
of T ~ 0.0004 from roughly a quarter of the BZ (100 
ky points) centred on k|| = Y. The total transmission, 
Ttotai = T s +T d = 0.93 + 0.04 = 0.97, compared to a 
transmission of 0.99 in the absence of disorder. Similar 
results were obtained for other k|| points. In the major- 
ity case, there is thus a strong specular peak surrounded 
by a weak diffuse background. 

The minority-spin Fermi surface projections of fee Cu 
and Co are shown in Figs. Illf al and lllf bh respectively. 
Compared to the corresponding panels in Fig. the 20 x 
20 point representation is seen to be sufficient to resolve 
the individual Fermi surface sheets of Co. To study the 
effect of interface disorder, we consider scattering out of 
two different k||S in Cu (Tigs. HTT c) and (d)). The first 
thing to note is the similarity of both transmission plots 
to the projected FS of Co, Fig. ITTT bl. suggesting very 
strong diffusive scattering proportional to the density of 
available final states. 

The first case we consider is where ky — Y for which 
the transmission was zero as a result of the symmetry 
of the states along the k y axis in the absence of disor- 
der. T(k|| = y, kj| ) is shown in Fig. llOf c). By con- 
trast with the majority-spin case just examined, there 
is now scattering to all other k-points in the 2D BZ, 




FIG. 11: Fermi surface projections of minority-spin fee Cu (a) 
and Co (b) derived from a single k-point using a 20 x 20 lateral 
supercell. The dark red point in the Cu Fermi surface projec- 
tion corresponds to the point Y' in the top righthand panel 
of Fig. © (c) T(F,k'|) and (d) T(F',k'|) calculated using 
20 different disorder configurations; the ballistic component 
T(Y' , ky = Y') is indicated by a white point because its value 
goes off scale. Results are for a "standard" configuration, 59 
20 different configurations of disorder were used. 
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Ek'^kJ^W = Y,-k\ /7) = 0.58 while T(Y, Y) has 
only increased from 0.00 in the clean case, to 0.01 in the 
presence of disorder. The effect of disorder is to increase 
the total transmission, T tota i(Y) = J2k' ^(^11 = ^'Kl^ 
from 0.00 to T S (Y) + T d (Y) = 0.01 + 0.58 = 0.59; for 
states which were originally strongly reflected, disorder 
increases the transmission. 

The second case we consider is that of a k-point slightly 
further away from the origin A along the k y axis which 
had a high transmission, T(Y') = 0.98, in the absence 
of disorder. For this k-point, T(k|| = F',k||), shown 

in Fig. HOf d). looks very similar to Fig. HOf ch There 
is strong diffuse scattering with J2k' ^k = Y', k'| ^ 

Y') = 0.54 while T(Y', Y') has been drastically decreased 
from 0.98 in the clean case, to 0.06 as a result of disorder. 
The total transmission, T total {Y') = T S {Y') + T d {Y') = 
0.06 + 0.54 = 0.60, is almost identical to what was found 
for the Y point. The effect of disorder has been to de- 
crease the transmission for states which were originally 
weakly reflected. The strong k-dependence of the trans- 
mission found in the specular case is largely destroyed by 
a small amount of disorder in the minority-spin channel. 
The contribution from specular component (integrated 
over 2D BZ) is reduced to 15% of the total transmission. 



E. Interface resistance 

To the best of our knowledge, spin-dependent interface 
transmissions have not yet been measured directly. What 
is usually don o 70 i 71 is to measure total resistances for a 
whole series of magnetic multilayers in which the total 
number of interfaces and/or the thicknesses of the indi- 
vidual layers is varied. The measured results are inter- 
preted in terms of volume resistivities and interface resis- 
tances. By applying an external magnetic field, the mag- 
netizations of neighbouring layers which are oriented an- 
tiparallel (AP) can be forced to line up in parallel (P). By 
measuring the resistances in both cases, spin-dependent 
volume resistivities and interface resistances can be ex- 
tracted using the two current series resistor modeliS2^2iS£ 
If we take expression (|48|l which relates the interface 
transmission to the interface resistance occurring in the 
2CSR model as given^^ we can study how typical uncer- 
tainties in interface transmission, arising from arbitrary 
assumptions about the interface disorder, lattice constant 
or basis set translate into uncertainty in predicted in- 
terface resistances. Using the transmission probabilities 
from Fig.Elin (|48J) results in the curves shown in Fig. 1121 
For comparison, a range of literature values for the spin- 
dependent interface resistances derived from experiments 
on sputtered and MBE (molecular beam exitaxy) grown 
multilayers 7 ^ is included in the figure. 

For the minority-spin case, experimental values 
(in units of film 2 ) range from 1.30-1.80 compared 
to calculated values of 1.29 for Cu[Cu.3Co.7]Co, 
through 1.37 for a disorder-free interface, to a 
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FIG. 12: Interface resistance for disordered interfaces as a 
function of the alloy concentration used to model disordered 
interfaces calculated using IJ48H and the transmission proba- 
bilities shown in Fig. 03 The experimental values for sputtered 
and MBE grown multilayers cited in Table I of Ref. |23 span 
a range of values which is indicated by the shaded regions. 



value of 2.25 for the 4ML model with x = 0.5, 

Cu[Cu.83Co.i7|Cu.67Co.33|Cu.33Co.67|Cu.i7Co. 83 ]Co. 

The influence of lattice constant and basis set on the 
clean interface resistance values is small (see Table . 
The present modelling of interface alloying shows that 
the interface resistance is more strongly dependent on 
the detailed spatial distribution of disorder than was 
previously found— where only the concentration range 
x = 0.5 ±0.06 of the 2ML interface alloy model extracted 
from experimenljSSiSiSiSi was explored. 

For the majority-spin case, the spread in values of 
the interface resistance extracted from experiment (for 
the same samples as for the minority-spin case) is quite 
small, 0.22-0.25, and does not overlap with the values 
of 0.34 found for a lattice constant of a = 3.614A Un- 
like the minority-spin case, changing the lattice constant 
or using an spdf basis leads to substantially larger val- 
ues (Table ^Qj. Because the majority-spin transmission 



a(A) 


3.549 




3.614 


Basis 


spdf 


spd 


spd 


i? maJ (lll) 


0.46 


0.39 


0.34 


i? min (lll) 


1.33 


1.32 


1.37 



TABLE III: Interface resistances, in units of ff2m 2 , for or- 
dered interfaces, calculated using expression 148H and the data 
from Tables [I] and [H] The values given here for a lattice 
constant of a = 3.614A differ slightly from those reported 
in Ref. |3 which were performed using energy-independent 
muffin-tin orbitals linearized about the centers of gravity of 
the occupied conduction states and not at the Fermi energy. 
The current implementation-*- uses energy-dependent, (non- 
linearized) MTO's, calculated exactly at the Fermi energy 
which improves the accuracy at no additional cost. 
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FIG. 13: Differential interface resistance as the number of in- 
terfaces increase for a disordered Cu/Co(lll) multilayer em- 
bedded between Cu leads. A 10 x 10 lateral supercell was 
used and the interface was modelled as two layers of 50%- 
50% alloy (2ML model). The results represent an average 
over 5 disorder configurations and were obtained for a "stan- 
dard" configuration,*— The range of experimental values— is 
indicated by the shaded regions. 



does not depend on the details of the interface disor- 
der, this cannot be the origin of the discrepancy. Mo- 
tivated by the weak scattering in this case, we examine 
the validity 69 i 72 i 73 i 74 i 75 of the 2CSR model by calculating 
the resistance of a magnetic multilayer containing a large 
number of disordered interfaces and plot the resistance 
added by each additional interface in Fig. ED Compared 
to similar calculations in Ref.0, the number of interfaces, 
size of lateral supercell (10 x 10) and disorder configura- 
tions averaged over are increased substantially. While the 
calculations are in very good agreement with Ohm's law 
for the strongly scattering minority-spin case, it can be 
seen that this is not the case for the majority-spin elec- 
trons. For a small number of interfaces there is a clear 
breakdown of Ohm's law and thus of the 2CSR model. 
The interface resistance eventually saturates at a value 
much lower than those extracted from experiment. While 
inclusion of bulk scattering will modify this picture some- 
what, exploratory calculations 76 indicate that the type of 
"bulk" impurities which may be reasonably expected to 
be found in sputtered or MBE grown multilayers affect 
the minority spin electrons much more than the majority 
spins. Agreement for the latter can only be achieved at 
the expense of ruining good agreement for the former. 



IV. DISCUSSION 

Details of a muffin tin orbital-based method suit- 
able for calculating from first-principles scattering ma- 
trices involving layered magnetic materials have been 
given. In a wide range of applicationSfSi^^LiSi^iL 5 ^ it 
has been shown to be much more efficient and transpar- 



ent than a previously used LAPW-based methodA 6 ^ 
Various other schemes have been developed for calcu- 
lating the transmission of electrons through an interface 
(or a more extended scattering region) both from first 
principles, 7 ! 8 ! 10 ! 11 ! 12 ! 14 ! 15 !! 7 ! 18 ! 35 ! 36 ! 38 ! 77 or using as input 
electronic structures which were calculated from first 
principles— i* 2 ^ 2 ^* 7 ^ 7 ^ ^ 8 ^ Most are based upon a for- 
mulation for the conductance in terms of non-equilibrium 
Green's functions^ 2 (NEGF) which reduces in the appro- 
priate limit to the well known Fisher-Lee (FL) linear- 
response formS 7 . for the conductance of a finite disor- 
dered wire embedded between crystalline leads. Most 
implementations of the NEGF or FL schemes have two 
disadvantages, (i) The transmission is calculated for a 
complex energy which leads to difficulties in studying 
for example, tunneling magnetoresistance, where the fi- 
nite imaginary part can give rise to an exponential de- 
cay which obscures the interesting physical decay of the 
transmission as a function of the barrier thickness, (ii) 
For a given value of transverse crystal momentum, the 
transmission is expressed as a trace over the basis set in 
terms of which the Green's function and self-energy are 
expressed. 54 While this has the advantage that the total 
transmission can be calculated without explicitly deter- 
mining the scattering states and can be computationally 
efficient, summation of the contributions from multiple 
scattering states can obscure real physical effects, for ex- 
ample, the role of the symmetries of individual scattering 
states seen in Fig. Explicit determination of the scat- 
tering states not only makes a detailed analysis of the 
scattering possible. The full scattering matrix, expressed 
in terms of the scattering states, can be used to bridged 
the gap between first-principles electronic structure cal- 
culations and phenomcnological models of transport used 
to analyse complex situations where a full first-principles 
treatment is not practical. 

We have instead made use of an alternative technique, 
suitable for Hamiltonians that can be represented in 
tight-binding form, that was formulated by Ando* 5 - and is 
based upon direct matching of the scattering-region wave 
function to the Bloch modes of the leads. The relation- 
ship between the wave function matching 4 ^ and Green 
function^LS approaches is not immediately obvious. It 
was suggested recently that WFM was incomplete^ 3 , but 
the equivalence of the two approaches could be proven^ 
Schemes similar in spirit to our own, but based upon em- 
pirical tight-binding Hamiltonians have been presented 
by Sanvito et al*2> and by Velev— i2i In contrast to 
these schemes, our TB-MTO formalism is a parameter- 
free approach that has all of the advantages derived from 
self-consistent determination of potentials and spin den- 
sities for systems for which these are not known from 
experiment. Judging from the size of systems to which 
it has been applied, it would seem that our implemen- 
tation is nevertheless substantially more efficient than 
these empirical schemes. The scattering regions treated 
in Figs. I7I9I10I and ITT1 contained as many as 3200 atoms 
(20 x 20 lateral supercell x 8 principal layers where the 



19 



potential was allowed to deviate from its bulk values) 
or, in the case of Fig. U31 ~ 15000 atoms (10 x 10 
lateral supercell x 150 principal layers). Our WFM 
scheme should not be confused^ with a recently devel- 
oped transport formalim 8 *^ also based upon TB-LMTOs 
but which makes use of the Caroli NEFG expression for 
the conductance in terms of a trace and a complex en- 
ergy. Khomyakov and Brocks^ have developed a scheme 
analogous to ours but based upon pseudopotentials and a 
real space grid which make it more suitable for studying 
quantum wires or the type of open structures studied in 
molecular electronics, but is computationally much more 
expensive. 

A third approach based upon "embedding"^ 4 *^ has 
been combined with full-potential linearized augmented 
plane wave method to yield what is probably the most 
accurate scheme to dat o 14 ' 1 ^?? but like the real space grid 
WFM method^ these methods are numerically very de- 
manding. 



APPENDIX A: VELOCITIES 

Expressions for the velocities of the propagating modes 
in the leads (Sect. IiTB|| are more easily derived using an 
energy independent Hamiltonian than the energy depen- 
dent tail-cancellation condition of section Hi Al To do so, 
we make use of the close relationship between the KKR 
tail-cancellation equation (|15fl and the linearized MTO 
(LMTO) Hamiltonian, both of which can be expressed 
in terms of the Hermitian matrisaia^ 

h a (e) = -[P Q (e)]- 1/2 (P Q (e) - S a ) [P a {e)]- 1/2 
= -P a (e)[P a {e)}- 1 + [P a {e)]- l/2 S a [P a {e)]- 1/2 . (Al) 

Fixing the energy at e = e F and defining the potential 
parameters" 

y/d^ = [P a (s F )}~ 1/2 (A2a) 



V. SUMMARY 

Details of a wave-function matching method suitable 
for calculating the scattering matrices in magnetic metal- 
lic hybrid structures based upon first-principles tight- 
binding muffin tin orbitals have been given and illus- 
trated with calculations for a variety of Co/Cu(lll) 
interface-related problems. The minimal basis of local- 
ized orbitals is very efficient, allowing large lateral su- 
percells to be handled. This allow us to model materials 
with large lattice mismatch or to study transport in the 
diffusive regime. Because the scattering states are calcu- 
lated explicitly, the effect of various types of scattering 
can be analyzed in detail. 
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c a = -P a {e F )/P a {e F )+e F , (A2b) 

(IA1|) can be written as 

h a = h a {e F )=c a + Vd" S a Va¥~e F (A3) 

Equation l|A3|) has the form of a two-center tight bind- 
ing Hamiltonian whose energy is given relative to e F . 
It provides the lowest order approximation 41 *' to the 
full LMTO Hamiltonian and yields eigenvalues correct 
to first order in (e — e F ). For eigenvalues equal to e F , it 
yields eigenvectors which are equal to those determined 
by the tail-cancellation condition (|15|) . up to a scaling 
factor (F")" 1 / 2 . 

To calculate the group velocities of states precisely 
at the linearization energy, in the present case at the 
Fermi energy, the first-order Hamiltonian (|A3|) can be 
used since any error vanishes identically for e(k) = e F . 
Using the translational symmetry of the leads, the Hamil- 
tonian (|A3|) for Bloch vector k is 

h^B/v W = £ e ' kTh RLAR' +TW ( A4 ) 

T 

where RL labels the sites and orbitals within the unit cell 
and T runs over lattice vectors. The energy eigenvalues 
e M (k) are the expectation values 

s M (k) = at (k)A a (k)a„(k) (A5) 

where the eigenvectors a M (k) are indexed by RL and we 
assumed normalization a^ • a^ = 1. It is now straightfor- 
ward to calculate the group velocity of the propagating 
mode 

/ , a RL^ l RL,(R'+T)L' a R'L l 
RL.R'L' 
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In the mixed representation I/, ky^ denned in section lil Bl 
(|A6() gives for the velocity in the stacking direction 

v» = j [a^?, /+ i(k||)A^-h.c] (A7) 

where d is the distance between equivalent monolayers in 
adjacent principal layers (PL), the hopping is assumed 
(as in section Til Bl) to extend only between neighbouring 
PLs and A M = exp(ik ■ T°) with T° connecting equiva- 
lent sites in the neighbouring PLs. Using the definition 
(|A1|) of h a and recalling that the solutions u M of the tail- 
cancellation equation (|15|l take implicitly into account 
the scaling factor (P a ) _1 / 2 we arrive at equation (|37() . 

APPENDIX B: SYMMETRY RELATIONS 

If we look closely at the transmission probabilities in 
Fig. El we see that the sheet resolved transmissions ex- 
hibit the geometrical symmetry of the underlying lattice 
(i.e. the three-fold rotational axis). The total trans- 
mission probability on the other hand possesses an extra 
inversion symmetry, T(k||) = T(— k||), which results in 
plots with a six- fold rotational axis. This higher symme- 
try is the manifestation of the fundamental time-reversal 
symmetry obeyed in the absence of spin-orbit coupling 
and a magnetic field. In the case of the bulk system 
time-reversal symmetry grants that for every eigenstate 
ip a (k) there exists the counterpart with the same energy 
and opposite wave vector (i.e. e Q (k) = s a (— k)) and 
the wave functions are related by the complex conjugate. 
The situation is more complicated in the case of the scat- 
tering state. Consider a state incoming from the left lead 
and scattered in the middle region. The wave function 
consists then of the incoming and reflected states in the 
left lead 

*£(k u ) = ^+(k||) + ^v AI (k||)v;,(k||) (bi) 

and of the transmitted states in the right lead 

^(k||)=^^(k||^+(k||). (B2) 

V 

The time reversal operation transforms the above "re- 
tarded" state into the "advanced" one in which a number 



of incoming states (from the left and the right) combine 
to produce a single outgoing state on the left, i.e. 

*£(-k,|) = ^r^(kn)V+(-k||) +^(-k,|) (B3) 
and 

^(-k||)=^^(k||)^-(-k||). (B4) 

V 

Equations (|B3(1 and 1|B4(1 impose a set of conditions on 
the values of scattering coefficients for the states with 
— k||. Combined with the analogous conditions derived 
for the states with the incoming state in the right lead, 
they are compactly expressed as 

7 = S(-k||)S*(k||) => 5(-k,|) =5 T (k,|). (B5) 

The scattering matrix S is defined as 

where A'^ and are matrices in the space of the lead 
modes and the primed coefficients describe scattering of 
the states incoming from the right. More specifically we 
have: 

W(- k ||) = *^( k ||) and VM(- k ||) = r w'( k ||) ( B7 ) 
Equation (|B7|I gives 

TM-h) = E IM- k n)l 2 - E M k n)l 2 = T ^( k n) 

(B8) 

In addition, for any two-terminal device, the Hermitic- 
ity of the scattering matrix guarantees that T-jic(k\\) = 
T£7j(k||) (see Ref. l26|) which finally proves the in-plane 
inversion symmetry mentioned at the beginning. The 
last step can not however be taken for the partial (FS re- 
solved) transmission probabilities. These quantities thus 
possess only the geometrical symmetry of the system. 
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